library(MASS)
library(Zelig)
library(xtable)
library(memisc)
source("C:/Users/Stephen/Documents/R/win-library/2.13/memisc/R/mtable-ext.R")

setwd("C:/WTO_Midwest/")

WTOdata <- read.csv(file="WTO_cox.csv",head=TRUE,sep=",")
attach(WTOdata)
names(WTOdata)

### Models WITHOUT Month Trends
### Base Model: PE1, UE, PE1UE, Retal
z.m1 <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m1

### Base Model + Plaintiff PCGDP
z.m2 <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m2

### Full Model: Base Model + Plaintiff PCGDP + Plaintiff Dynamics
z.m3 <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP1_PlUE + Plaintiff_UE + PlE1,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m3

### Full Model: Base Model + Plaintiff PCGDP + Plaintiff Dynamics + LegalSum
z.m4 <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP1_PlUE + Plaintiff_UE + PlE1 + LegalSum,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m4


### Models WITH Quadratic Month Trends
###
z.m1m <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m1m

### Base Model + Plaintiff PCGDP
z.m2m <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m2m

### Full Model: Base Model + Plaintiff PCGDP + Plaintiff Dynamics
z.m3m <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP1_PlUE + Plaintiff_UE + PlE1 + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m3m

### Full Model: Base Model + Plaintiff PCGDP + Plaintiff Dynamics + LegalSum
z.m4m <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP1_PlUE + Plaintiff_UE + PlE1 + LegalSum + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m4m



### Table of Cox Results
toLatex(
  relabel(
       mtable("Model 1"=z.m1,
              "Model 2"=z.m3,
              "Model 3"=z.m1m,
		  "Model 4"=z.m3m
              )
          ),
  ddigits=5)


### Models with Efron, Exact
z.m1EF <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("efron"))
z.m1EF
z.m3EF <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP1_PlUE + Plaintiff_UE + PlE1,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("efron"))
z.m3EF
z.m1mEF <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("efron"))
z.m1mEF
z.m3mEF <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP1_PlUE + Plaintiff_UE + PlE1 + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("efron"))
z.m3mEF


















### Full Model b: Base Model + Plaintiff PCGDP + Plaintiff Dynamics b
z.m3mb <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP_UE + Plaintiff_UE + Elec_Prox_Plaintiff + Month + Month2,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m3mb

### Full Model b: Base Model + Plaintiff PCGDP + Plaintiff Dynamics b + LegalSum
z.m4mb <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_UE + Unemp_6mo + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + PlEP_UE + Plaintiff_UE + Elec_Prox_Plaintiff + Month + Month2 + LegalSum,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"))
z.m4mb





###
# Survival Curves
###

###
# This sets up the data, for PE1/0 and UE3-7
se1pe1ue3<-WTOdata[0:1,]
se1pe1ue3$PE1<-as.matrix(rep(1,1), nrow=1, ncol=1)
se1pe1ue3$PE1_UE<-as.matrix(rep(3,1), nrow=1, ncol=1)
se1pe1ue3$Unemp_6mo<-as.matrix(rep(3,1), nrow=1, ncol=1)
se1pe1ue3$Exports_6mo<-as.matrix(rep(mean(WTOdata$Exports_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$Imports_6mo<-as.matrix(rep(mean(WTOdata$Imports_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$Plaintiff_PCGDP<-as.matrix(rep(mean(WTOdata$Plaintiff_PCGDP,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$Plaintiff_UE<-as.matrix(rep(mean(WTOdata$Plaintiff_UE,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$PlEP1_PlUE<-as.matrix(rep(mean(WTOdata$PlEP_UE,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$PlE1<-as.matrix(rep(1,1), nrow=1, ncol=1)
se1pe1ue3$Month<-as.matrix(rep(mean(WTOdata$Month,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$Month2<-as.matrix(rep(mean(WTOdata$Month2,na.rm=TRUE),1), nrow=1, ncol=1)
#se1pe1ue3$Month3<-as.matrix(rep(mean(WTOdata$Month3,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$PlEP_UE<-as.matrix(rep(mean(WTOdata$PlEP_UE,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$Elec_Prox_Plaintiff<-as.matrix(rep(mean(WTOdata$Elec_Prox_Plaintiff,na.rm=TRUE),1), nrow=1, ncol=1)
se1pe1ue3$LegalSum<-as.matrix(rep(mean(WTOdata$LegalSum,na.rm=TRUE),1), nrow=1, ncol=1)



se1pe1ue4<-se1pe1ue3
se1pe1ue4$PE1_UE<-as.matrix(rep(4,1), nrow=1, ncol=1)
se1pe1ue4$Unemp_6mo<-as.matrix(rep(4,1), nrow=1, ncol=1)

se1pe1ue5<-se1pe1ue3
se1pe1ue5$PE1_UE<-as.matrix(rep(5,1), nrow=1, ncol=1)
se1pe1ue5$Unemp_6mo<-as.matrix(rep(5,1), nrow=1, ncol=1)

se1pe1ue6<-se1pe1ue3
se1pe1ue6$PE1_UE<-as.matrix(rep(5,1), nrow=1, ncol=1)
se1pe1ue6$Unemp_6mo<-as.matrix(rep(5,1), nrow=1, ncol=1)

se1pe1ue7<-se1pe1ue3
se1pe1ue7$PE1_UE<-as.matrix(rep(5,1), nrow=1, ncol=1)
se1pe1ue7$Unemp_6mo<-as.matrix(rep(5,1), nrow=1, ncol=1)

se1pe0ue3<-se1pe1ue3
se1pe0ue3$PE1_UE<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1pe0ue3$PE1<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1pe0ue4<-se1pe0ue3
se1pe0ue4$Unemp_6mo<-as.matrix(rep(4,1), nrow=1, ncol=1)
se1pe0ue5<-se1pe0ue3
se1pe0ue5$Unemp_6mo<-as.matrix(rep(5,1), nrow=1, ncol=1)
se1pe0ue6<-se1pe0ue3
se1pe0ue6$Unemp_6mo<-as.matrix(rep(6,1), nrow=1, ncol=1)
se1pe0ue7<-se1pe0ue3
se1pe0ue7$Unemp_6mo<-as.matrix(rep(7,1), nrow=1, ncol=1)


pe1ue3<-as.data.frame(se1pe1ue3[1,])
pe1ue4<-as.data.frame(se1pe1ue4[1,])
pe1ue5<-as.data.frame(se1pe1ue5[1,])
pe1ue6<-as.data.frame(se1pe1ue6[1,])
pe1ue7<-as.data.frame(se1pe1ue7[1,])
pe0ue3<-as.data.frame(se1pe0ue3[1,])
pe0ue4<-as.data.frame(se1pe0ue4[1,])
pe0ue5<-as.data.frame(se1pe0ue5[1,])
pe0ue6<-as.data.frame(se1pe0ue6[1,])
pe0ue7<-as.data.frame(se1pe0ue7[1,])



### This generates the survival curve figure, using Model 3 with trends
# PE1 and UE4/6
par(mfrow=c(2,1),oma=c(.3,.3,2,.3))
plot(survfit(z.m4m, newdata=pe1ue6), xlim=c(0,100), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Unemployment = 6%")
plot(survfit(z.m4m, newdata=pe1ue4), xlim=c(0,100), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Unemployment = 4%")
title("Survival Curves for WTO Dispute in Election Year",
	 outer=TRUE, cex=.3)
dev.print(postscript, "C:/WTO_Midwest/Drafts/PE1UE4_6.eps", horizontal=FALSE)









par(mfrow=c(2,2))
plot(survfit(z.out.m1, newdata=pe1ue4), ylab="survival proportion", xlab="month (age)", main="PE1, UE4")
plot(survfit(z.out.m1, newdata=pe1ue6), ylab="survival proportion", xlab="month (age)", main="PE1, UE6")
plot(survfit(z.out.m1, newdata=pe0ue4), ylab="survival proportion", xlab="month (age)", main="PE0, UE4")
plot(survfit(z.out.m1, newdata=pe0ue6), ylab="survival proportion", xlab="month (age)", main="PE0, UE6")


plot(survfit(z.out.m1, newdata=pe1ue3), ylab="survival proportion", xlab="month (age)", main="PE1, UE3")
plot(survfit(z.out.m1, newdata=pe1ue5), ylab="survival proportion", xlab="month (age)", main="PE1, UE5")
plot(survfit(z.out.m1, newdata=pe0ue3), ylab="survival proportion", xlab="month (age)", main="PE0, UE3")
plot(survfit(z.out.m1, newdata=pe0ue4), ylab="survival proportion", xlab="month (age)", main="PE0, UE4")
plot(survfit(z.out.m1, newdata=pe0ue5), ylab="survival proportion", xlab="month (age)", main="PE0, UE5")
plot(survfit(z.out.m1, newdata=pe1ue6), ylab="survival proportion", xlab="month (age)", main="PE1, UE6")
plot(survfit(z.out.m1, newdata=pe0ue6), ylab="survival proportion", xlab="month (age)", main="PE0, UE6")
plot(survfit(z.out.m1, newdata=pe0ue7), ylab="survival proportion", xlab="month (age)", main="PE0, UE7")


# These are the two figures currently in the JMS for PE1/0 and UE4/6
par(mfrow=c(2,1),oma=c(.3,.3,2,.3))
plot(survfit(z.out.m1, newdata=pe1ue4), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Unemployment = 4%")
plot(survfit(z.out.m1, newdata=pe1ue6), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Unemployment = 6%")
title("Survival Curves for WTO Dispute in Election Year",
	 outer=TRUE, cex=.3)
dev.print(postscript, "C:/Users/Stephen/Desktop/JMP/PE1UE4_6.eps", horizontal=FALSE)

par(mfrow=c(2,1),oma=c(.3,.3,2,.3))
plot(survfit(z.out.m1, newdata=pe0ue4), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Unemployment = 4%")
plot(survfit(z.out.m1, newdata=pe0ue6), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Unemployment = 6%")
title("Survival Curves for WTO Dispute in Non-election Year",
	 outer=TRUE, cex=.3)
dev.print(postscript, "C:/Users/Stephen/Desktop/JMP/PE0UE4_6.eps", horizontal=FALSE)

###
# Substantive Effects of Exports
# Vary EX from .008 (25th percentile) to .05 (75th percentile), PE=1
se1ex008<-se1pe1ue3
se1ex008$Unemp_6mo<-as.matrix(rep(mean(WTOdata$Unemp_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1ex008$PE1_UE<-as.matrix(rep(mean(WTOdata$Unemp_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1ex008$Exports_6mo<-as.matrix(rep(0.008,1), nrow=1, ncol=1)
se1ex05<-se1ex008
se1ex05$Exports_6mo<-as.matrix(rep(0.05,1), nrow=1, ncol=1)
ex008<-as.data.frame(se1ex008[1,])
ex05<-as.data.frame(se1ex05[1,])
par(mfrow=c(2,1),oma=c(.3,.3,2,.3))
plot(survfit(z.out.m1, newdata=ex008), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Exports = 0.08%")
plot(survfit(z.out.m1, newdata=ex05), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Exports = 5%")
title("Survival Curves for WTO Dispute",
	 outer=TRUE, cex=.3)
dev.print(postscript, "C:/Users/Stephen/Desktop/JMP/EX008_05.eps", horizontal=FALSE)

###
# Substantive Effects of Imports
# Vary IM from .009 (25th percentile) to .075 (75th percentile), PE=1
se1im009<-se1pe1ue3
se1im009$Unemp_6mo<-as.matrix(rep(mean(WTOdata$Unemp_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1im009$PE1_UE<-as.matrix(rep(mean(WTOdata$Unemp_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1im009$Imports_6mo<-as.matrix(rep(0.009,1), nrow=1, ncol=1)
se1im075<-se1im009
se1im075$Imports_6mo<-as.matrix(rep(0.075,1), nrow=1, ncol=1)
im009<-as.data.frame(se1im009[1,])
im075<-as.data.frame(se1im075[1,])
par(mfrow=c(2,1),oma=c(.3,.3,2,.3))
plot(survfit(z.out.m1, newdata=im009), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Imports = 0.09%")
plot(survfit(z.out.m1, newdata=im075), ylab="survival proportion", xlab="age of tariff (in months)", main="U.S. Imports = 7.5%")
title("Survival Curves for WTO Dispute",
	 outer=TRUE, cex=.3)
dev.print(postscript, "C:/Users/Stephen/Desktop/JMP/IM009_075.eps", horizontal=FALSE)


###
# Substantive Effects of Plaintiff PCGDP
# Vary EX from .008 (25th percentile) to .05 (75th percentile), PE=1
se1p25<-se1pe1ue3
se1p25$Unemp_6mo<-as.matrix(rep(mean(WTOdata$Unemp_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1p25$PE1_UE<-as.matrix(rep(mean(WTOdata$Unemp_6mo,na.rm=TRUE),1), nrow=1, ncol=1)
se1p25$Plaintiff_PCGPD<-as.matrix(rep(4000,1), nrow=1, ncol=1)
se1p75<-se1p25
se1p75$Plaintiff_PCGPD<-as.matrix(rep(27000,1), nrow=1, ncol=1)
p25<-as.data.frame(se1p25[1,])
p75<-as.data.frame(se1p75[1,])
par(mfrow=c(2,1),oma=c(.3,.3,2,.3))
plot(survfit(z.out.m1, newdata=p25), ylab="survival proportion", xlab="age of tariff (in months)", main="Plaintiff Per Capita GDP = $4,000")
plot(survfit(z.out.m1, newdata=p75), ylab="survival proportion", xlab="age of tariff (in months)", main="Plaintiff Per Capita GDP = $27,000")
title("Survival Curves for WTO Dispute",
	 outer=TRUE, cex=.3)
dev.print(postscript, "C:/Users/Stephen/Desktop/JMP/PCGDP25_75.eps", horizontal=FALSE)





#mtext("Survival curves holding other variables at sample means, using results from Model 3.",
#	 SOUTH<-1, line=0.1, adj=0, cex=.6, outer=TRUE)
#mtext("Dotted lines are 95% confidence intervals.",
#	 SOUTH<-1, line=0.6, adj=0, cex=.6, outer=TRUE)













# Coxph Model 2: PE1 + TPO_imp
z.out.m2 <- zelig(Surv(age_start,age_stop,fstatus3) ~ PE1_TPOimp + TPO_imp + PE1 + Exports_6mo + Imports_6mo + Plaintiff_PCGDP + Plaintiff_UE + PlEP_UE + Elec_Prox_Plaintiff + LegalSum,
	model = "coxph", data = WTOdata, robust=TRUE, method=c("breslow"), control=t)
z.out.m2


se2pe1tpo46<-se1pe1ue3
se2pe1tpo46$TPO_imp<-as.matrix(rep(46,1), nrow=1, ncol=1)
se2pe1tpo46$PE1_TPOimp<-as.matrix(rep(46,1), nrow=1, ncol=1)
se2pe1tpo48<-se2pe1tpo46
se2pe1tpo48$TPO_imp<-as.matrix(rep(48,1), nrow=1, ncol=1)
se2pe1tpo48$PE1_TPOimp<-as.matrix(rep(48,1), nrow=1, ncol=1)
se2pe1tpo50<-se2pe1tpo46
se2pe1tpo50$TPO_imp<-as.matrix(rep(50,1), nrow=1, ncol=1)
se2pe1tpo50$PE1_TPOimp<-as.matrix(rep(50,1), nrow=1, ncol=1)
se2pe1tpo52<-se2pe1tpo46
se2pe1tpo52$TPO_imp<-as.matrix(rep(52,1), nrow=1, ncol=1)
se2pe1tpo52$PE1_TPOimp<-as.matrix(rep(52,1), nrow=1, ncol=1)

se2pe0tpo46<-se2pe1tpo46
se2pe0tpo46$PE1<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo46$PE1_TPOimp<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo48<-se2pe1tpo48
se2pe0tpo48$PE1<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo48$PE1_TPOimp<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo50<-se2pe1tpo50
se2pe0tpo50$PE1<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo50$PE1_TPOimp<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo52<-se2pe1tpo52
se2pe0tpo52$PE1<-as.matrix(rep(0,1), nrow=1, ncol=1)
se2pe0tpo52$PE1_TPOimp<-as.matrix(rep(0,1), nrow=1, ncol=1)

pe1tpo46<-as.data.frame(se2pe1tpo46[1,])
pe1tpo48<-as.data.frame(se2pe1tpo48[1,])
pe1tpo50<-as.data.frame(se2pe1tpo50[1,])
pe1tpo52<-as.data.frame(se2pe1tpo52[1,])
pe0tpo46<-as.data.frame(se2pe0tpo46[1,])
pe0tpo48<-as.data.frame(se2pe0tpo48[1,])
pe0tpo50<-as.data.frame(se2pe0tpo50[1,])
pe0tpo52<-as.data.frame(se2pe0tpo52[1,])

par(mfrow=c(2,2))
plot(survfit(z.out.m2, newdata=pe1tpo46), ylab="survival proportion", xlab="month (age)", main="PE1, TPO46")
plot(survfit(z.out.m2, newdata=pe1tpo52), ylab="survival proportion", xlab="month (age)", main="PE1, TPO52")
plot(survfit(z.out.m2, newdata=pe0tpo46), ylab="survival proportion", xlab="month (age)", main="PE0, TPO46")
plot(survfit(z.out.m2, newdata=pe0tpo52), ylab="survival proportion", xlab="month (age)", main="PE0, TPO52")



















par(mfrow=c(2,2),oma=c(2,1,2,1))
plot(survfit(z.out.CCm3, newdata=se3dembd0), ylab="survival proportion", xlab="quarter", main="Democracy, 0 Deaths")
plot(survfit(z.out.CCm3, newdata=se3ndembd0), ylab="survival proportion", xlab="quarter", main="Non-democracy, 0 Deaths")
plot(survfit(z.out.CCm3, newdata=se3dembd2k), ylab="survival proportion", xlab="quarter", main="Democracy, 2,000 Deaths")
plot(survfit(z.out.CCm3, newdata=se3ndembd2k), ylab="survival proportion", xlab="quarter", main="Non-democracy, 2,000 Deaths")
title("Survival Curves for ICC Ratification, by Regime Type and Conflict History",
	 outer=TRUE, cex=.6)
mtext("Survival curves holding other variables at sample means, using results from Model 3.",
	 SOUTH<-1, line=0.1, adj=0, cex=.6, outer=TRUE)
mtext("Dotted lines are 95% confidence intervals.",
	 SOUTH<-1, line=0.6, adj=0, cex=.6, outer=TRUE)
dev.print(postscript, "C:/Users/Stephen/Dropbox/ICC/Drafts/survcurves.eps", horizontal=FALSE)

par(mfrow=c(1,2))
plot(survfit(z.out.CCm3, newdata=se3dembd0), xlab="quarter", main="Democracy, 0 Deaths")
plot(survfit(z.out.CCm3, newdata=se3ndembd0), xlab="quarter", main="Non-democracy, 0 Deaths")
#plot(survfit(z.out.CCm3, newdata=se3dembd2k), xlab="quarter", main="Democracy, 2,000 Deaths")
#plot(survfit(z.out.CCm3, newdata=se3ndembd2k), xlab="quarter", main="Non-democracy, 2,000 Deaths")
#dev.print(eps, filename=("C:/Users/Stephen/Dropbox/ICC/Drafts/survcurves.eps"), width=1000, height=1000)



###
# This makes the survival curves for all the SD combinations
###
se1dem0<-ICCdata[0:1,]
se1dem0$cw_90<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1dem0$dem_med_hi_cw90<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1dem0$dem_med_hi<-as.matrix(rep(1,1), nrow=1, ncol=1)
se1dem0$logmilper1<-as.matrix(rep(mean(ICCdata$logmilper1,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$logpk<-as.matrix(rep(mean(ICCdata$logpk,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$corrected_regrat100<-as.matrix(rep(mean(ICCdata$corrected_regrat100,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$humrightstreat<-as.matrix(rep(mean(ICCdata$humrightstreat,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$extrconflict_ongoing<-as.matrix(rep(mean(ICCdata$extrconflict_ongoing,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$const__amend__required_for_ratif<-as.matrix(rep(mean(ICCdata$const__amend__required_for_ratif,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$icc_elected_officials<-as.matrix(rep(mean(ICCdata$icc_elected_officials,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$ldrs<-as.matrix(rep(mean(ICCdata$ldrs,na.rm=TRUE),1), nrow=1, ncol=1)
se1dem0$brit__leg__herit_1<-as.matrix(rep(mean(ICCdata$brit__leg__herit_1,na.rm=TRUE),1), nrow=1, ncol=1)

se1dem1<-se1dem0
se1dem1$cw_90<-as.matrix(rep(1,1), nrow=1, ncol=1)
se1dem1$dem_med_hi_cw90<-as.matrix(rep(1,1), nrow=1, ncol=1)

se1dem2<-se1dem0
se1dem2$cw_90<-as.matrix(rep(2,1), nrow=1, ncol=1)
se1dem2$dem_med_hi_cw90<-as.matrix(rep(2,1), nrow=1, ncol=1)

se1ndem0<-se1dem0
se1ndem0$dem_med_hi_cw90<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1ndem0$dem_med_hi<-as.matrix(rep(0,1), nrow=1, ncol=1)

se1ndem1<-se1dem1
se1ndem1$dem_med_hi_cw90<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1ndem1$dem_med_hi<-as.matrix(rep(0,1), nrow=1, ncol=1)

se1ndem2<-se1dem2
se1ndem2$dem_med_hi_cw90<-as.matrix(rep(0,1), nrow=1, ncol=1)
se1ndem2$dem_med_hi<-as.matrix(rep(0,1), nrow=1, ncol=1)

se1dem0<-as.data.frame(se1dem0[1,])
se1dem1<-as.data.frame(se1dem1[1,])
se1dem2<-as.data.frame(se1dem2[1,])
se1ndem0<-as.data.frame(se1ndem0[1,])
se1ndem1<-as.data.frame(se1ndem1[1,])
se1ndem2<-as.data.frame(se1ndem2[1,])

par(mfrow=c(3,2))
plot(survfit(z.out.SDm1, newdata=se1dem0), xlab="quarter", main="Democracy, CW 0")
plot(survfit(z.out.SDm1, newdata=se1ndem0), xlab="quarter", main="Nondemocracy, CW 0")
plot(survfit(z.out.SDm1, newdata=se1dem1), xlab="quarter", main="Democracy, CW 1")
plot(survfit(z.out.SDm1, newdata=se1ndem1), xlab="quarter", main="Nondemocracy, CW 1")
plot(survfit(z.out.SDm1, newdata=se1dem2), xlab="quarter", main="Democracy, CW 2")
plot(survfit(z.out.SDm1, newdata=se1ndem2), xlab="quarter", main="Nondemocracy, CW 2")
